Exploiting a semi-analytic approach to study first order phase transitions 
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In a previous contribution, Phys. Rev. Lett 107, 230601 (2011), we have proposed a method 
to treat first order phase transitions at low temperatures. It describes arbitrary order parameter 
through an analytical expression W , which depends on few coefficients. Such coefficients can be 
calculated by simulating relatively small systems, hence with a low computational cost. The method 
determines the precise location of coexistence lines and arbitrary response functions (from proper 
derivatives of W). Here we exploit and extend the approach, discussing a more general condition 
for its validity. We show that in fact it works beyond the low T limit, provided the first order phase 
transition is strong enough. Thus, W can be used even to study athermal problems, as exemplified 
for a hard-core lattice gas. We furthermore demonstrate that other relevant thermodynamic quan- 
tities, as entropy and energy, are also obtained from W. To clarify some important mathematical 
features of the method, we analyze in details an analytically solvable problem. Finally, we discuss 
different representative models, namely, Potts, Bell-Lavis, and associating gas-lattice, illustrating 
the procedure broad applicability. 
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I. INTRODUCTION 

First order phase transitions (FOPTs) occur in all sort 
of processes in nature [l[ , being extensively studied un- 
der different point of views and by a great diversity of 
approaches [2|. However, for strong FOPTs or at low 
temperatures, proper and reliable analysis still may be 
challenging [3|, |4(. This is so because in such contexts 
many simulation methods either can face technical diffi- 
culties (e.g., associated to very slow convergence) or de- 
mand considerable computational efforts (e.g., due to the 
necessity to simulate large systems). 

To overcome some of the above mentioned problems, 
in a recent short contribution [j| we have proposed a gen- 
eral semi-analytic method (a not so common approach in 
this area @) to deal with FOPTs at low T's. The method 
combines simple ideas, resulting in an accurate "combo" 
protocol to study FOPTs. Briefly (details in Sec. II): 
(a) considering a special decomposition for the partition 
function at low T's, an analytical expression W to char- 
acterize FOPTs (e.g., order parameter) is derived; (b) 
W depends on some coefficients, but which can be deter- 
mined through few numerical simulations (thus, from a 
computationally inexpensive procedure); and (c) highly 
profiting from the analyticity of W and using finite scale 
analysis for rather small systems, location of the tran- 
sition points, order parameters behavior, and response 
functions (like specific heats and compressibilities), are 
obtained with good precision. 

In the original paper [5|, we have demonstrated the 
framework power by means of different examples, includ- 
ing the analysis of complicated (and often hard to sim- 
ulate) Hamiltonians which describe diverse effects, such 
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as water-like anomalies and ferrimagnetic-ferromagnetic 
and ferromagnetic-ferromagnetic transitions. Neverthe- 
less, distinct important aspects of the approach were not 
addressed in-depth. 

Here we further explore the method, unveiling some of 
its mathematical and conceptual features. We address 
the approach extent of validity, demonstrating it can 
work fine beyond the originally derived regime of appli- 
cability, i.e., at low T's. We propose a concrete condition 
(testable by simple simulations) which shows it can lead 
to good results in instances of strong FOPTs. Besides 
order parameters and their derivatives (i.e., generalized 
susceptibilities), we discuss examples of other relevant 
thermodynamic quantities, like entropy and energy, that 
around the phase transition are also well described by 
W (given that in such cases, the coefficients for W are 
properly determined) . 

The paper is organized as the following. The method 
main ideas and key expressions are summarized in Sec. 
II (for completeness, full derivations are presented in the 
Appendix). To clarify important technical aspects of the 
approach, an exactly solvable model [7[ is discussed in 
Sec. III. In Sec. IV, an extensive analysis of repre- 
sentative models is carried out. The Potts model 
displaying extreme FOPTs for large q values (and for 
which the transition points are known exactly), is thor- 
oughly investigated, including entropy and energy. The 
method high numerical accuracy is illustrated with the 
Bell-Lavis model 0. Taking the associating lattice gas 
(ALG) model [10|,|ll| as an example, it is shown that the 
approach can be used for higher temperatures, provided 
the phase transition is sufficiently strong. Considering 
a hard-core gas lattice model, it is demonstrated that 
even athermal problems can be studied with the method. 
Guided by the numerical simulations and straightforward 
thermodynamic arguments, a general condition setting 
the approach validity is proposed in Sec. V. Numerically, 
it is based on the calculation of the order parameter mul- 
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timodal probability distribution at the coexistence. Fi- 
nally, remarks and conclusion are drawn in Sec. VI. 

II. THE METHOD 
A. Main ideas and results 

The method start point [5j is the fact that for finite 
systems at low temperatures and having J\f coexisting 
stable phases, the partition function is well described by 
a sum of J\f exponential terms, or [l2| (/3 = l/(kT)) 

A/" 

Z = J2 a n e M-PVf n ]. (1) 
n=l 

For each phase n, /„ is the free energy per volume V and 
a n the degeneracy (see also the Appendix). 

Typically, relevant thermodynamic quantities have the 
form W — — (f3V)~ 1 (d/dt;) ln[Z], with £ an appropriate 
control parameter (e.g., chemical potential /j,, tempera- 
ture, etc). For instance, if £ = /x, the density follows 
directly from p(/i,T) = ~W and if £ = fa the energy per 
volume is u = faW . 

As discussed in details in the Appendix, close to the 
transition point £ = £* and considering Eq. ([T]), one 
finds very generally that W can be approximated by (y = 

£-D 

w « (61 + h 

cxp[— Qij/])/(l + ^ a exp[— cijy]). (2) 

i=2 i=2 

The coefficients a,, bi and Ci are independent on the con- 
trol parameter £ and only the dj's are (linear) functions 
of V. In this way, at the coexistence (y = 0), W is inde- 
pendent on the volume and for all V the curves W x £ 
cross at £ = £*. Therefore, IF in Eq. ^ does not only 
describe order parameters, but it also gives the thermo- 
dynamic limit estimate for the transition point £*. 

Besides the above, two other aspects of the method, 
relevant for applications, are the following. First, the 
explicit dependence of Z on the free energy at low T's, 
Eq. ([1]). Close to £*, it makes both / and the entropy 
per volume s (once s = (u — f) /T) also to have the same 
functional form of Eq. Second, analytical derivatives 
of Eq. ©, e.g., leading to specific heat, susceptibility, 
and order parameters which are not necessarily first order 
derivatives of the free-energy, are easily calculated. 

Lastly, an important advantage of the present proce- 
dure is that Eq. ([3]) is valid for any volume. So, by 
considering relatively small V's we can obtain the pa- 
rameters a's, 6's, and c's with a low computational cost 
(see below). It allows to describe a first order phase tran- 
sition (at low temperatures) with a direct, accurate and 
numerically cheap method. Moreover, as we are also go- 
ing to show below, the approach in fact can be applied to 
broader situations than that initially assumed to derive 
Eq. ©, namely, of low T's. 



B. Numerical simulations 

Equation ([2]) is an analytical expression to describe 
proper order parameters (as well other thermodynamic 
functions) around the phase transition. Nevertheless, the 
coefficients a's, b's, and c's need to be determined. Al- 
though approximated expressions for these parameters do 
exist (Appendix) , much better results are obtained direct 
from numerics. The protocol is then: to use some sim- 
ulation method to generate the sought thermodynamic 
quantity for different values of £; to compare with the 
corresponding curve W; and to determine the coefficients 
by fitting. In particular, the actual general shape of W 
requires only few points for a proper adjustment, making 
the numerics rather fast. Finally, once the coefficients 
are known, finite size scale analysis, crossing determina- 
tion, calculation of derivatives, etc, can all be performed 
analytically. 

As a condition to choose any simulation approach to 
fit the coefficients of W, one should guarantee it is ap- 
propriate for the system at hands. Then, in this work we 
consider the parallel tempering (PT), which is general, 
simple to use and very efficient for low and intermediate 
system sizes, even for strong FOPTs [3, [HI, HH • Such fea- 
tures also qualify full simulations from the PT as good 
benchmarks to test Eq. 

For completeness, we briefly describe how to imple- 
ment the PT in our examples (for a very detailed dis- 
cussion, explaining each step in the method and its ap- 
plication to FOPT see, for instance, Ref. |3j). Basically, 
the PT (an enhanced sampling method) uses configura- 
tions from high to perform an ergodic walk at low T's. 
It simultaneously simulates a set of R replicas - in the 
temperature interval {T±, ...,Tr} - by means of a stan- 
dard algorithm (e.g., Metropolis, cluster, etc). When 
evolving any replica i at a temperature Tj (through an 
one- flip procedure), a given site k is chosen randomly 
and its state variable a' k may change to a new value a'l 
according to the probability = min{I, exp[— fa A"H]}, 
where AH = H(cr") — 7i(a') is the energy change due 
to the transition. Moreover, from time to time a pair of 
replicas (say, at and Tj and with microscopy configura- 
tions a' and a", respectively) can undergo a temperature 
switching, drawn from the probability 

P T ^ T] = min{l, exp[(A - ft) («(</) - H(a"))}}. (3) 

Typically, the number of replicas does not need to be 
very high. For the concrete calculations in this work we 
set R = 12. We also consider adjacent and non-adjacent 
replica exchanges. 

For the Potts model (Sec. IV. A), which presents strong 
discontinuous transition for large g's, we replace the 
above one-flip step by the Wolff cluster algorithm (l5j . 
In short, initially a seed site k (in the state <7fc) is chosen 
at will. Then, with probability p = (1 — exp[— j3J]) 5 ak ai 
(for J the two neighbor sites interaction energy when 
both are at o^) k is connected to each nearest neighbor 
site I. This is repeated for all the new sites of the cluster 
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until no extra site can be added. The entire cluster sites 
states are finally changed to a same a', randomly chosen 
from all the possible values for the state variable. 

For the athermal hard-core lattice model we shall ana- 
lyze here, very few modifications in the previous prescrip- 
tions are necessary. They are explained in Sec. IV. D. 

For all the examples in Section IV, we use only four 
points from the simulations to determine the coefficients 
in Eq. ([2]). The resulting analytical expressions are then 
compared with accurate numerics from the PT method 
described above. In many cases we also confront the 
present with different calculations in the literature to fur- 
ther check the efficiency of our general approach. 



III. AN EXACTLY SOLVABLE MODEL 

We begin with an analytical example, very instruc- 
tive to unveil certain mathematical aspects of the present 
method. Thus, we discuss a prototype model proposed 
in 0, which although simple, displays all the essential 
aspects associated to first order phase transitions. The 
problem grand-partition function is given by 0] 



Z=(l + z) v (l + z rV ) 



(4) 



where r is an arbitrary parameter and z = exp[/3 p] is the 
fugacity. In the thermodynamic limit of V — > oo, Z has 
a real root at z = 1 (Le., j3 p = 0), which according to 
the Yang-Lee theory [l6| characterizes a FOPT between 
two phases {N = 2) . 

Here, an appropriate order parameter is the density 



rz 



rV 



1 



1 + z 



rV ' 



(5) 



Indeed, for V — > oo p has a gap of r since in such limit 
p(z = 1~) = 1/2 and p(z = 1+) = 1/2 + r. Moreover, 
p(z = 1) = (1 + r)/2 regardless of V. 

Now, consider a finite p. So, low temperatures corre- 
spond to large values for z (large /3 p) and a good approx- 
imation for Eq. (jU) is Z sa exp[/3pV] + exp[/3(r + l)pV], 
which obviously is in the general form of Eq. (JT]). How- 
ever, the phase transition takes place at z = 1. For z rj 1 
in Eq. (QJ, we can consider 1 + z k, 2y/z, getting 



Zr*exp[-0Vf o ]+exp[-PVf r ], 



(6) 



with 



fa 



1 



((r + -)/3/i + ln[2]). (7) 



Notice that again Z takes the form of Eq. (JXJ) . Hence, 
although being a particular model, it shows that Eq. ([1]) 
and so Eq. @ can hold true in more general instances 
than just that of low T's (see the discussion in Sec. V). 

As mentioned in Sec. II, the density is given by p = 
—W with £ = p. Thus, using the analytical relations for 
the coefficients of Eq. © in the Appendix, one finds that 
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FIG. 1. p vs. /3p for distinct 7=1x1 and r = 1/2. The 
left inset is a blow up around the transition point, at which 
p = 3/4. Symbols correspond to the exact p and dashed lines 
to Eq. (JSJ). The continuous lines are for the density given by 
the general Eq. ([2]), with the parameters C2 and b's the same 
ones used in Eq. {SJ, but 02//? from a best numerical fitting. 
For two V values and r — 1/2, Po(p) and P r (p) are plotted 
together in the right inset. Since Pq (P r ) is practically zero for 
p > 3/4 (p < 3/4), the overlap between the two distributions 
is negligible. 



for the present case: c 2 = 1, b\ = —1/2, b 2 = — (r+ 1/2), 
and ai = r(3V. Therefore, close to the phase transition 
point it follows from Eq. ^ an approximation for p, or 



l + z 



rV ' 



(8) 



which yields the correct limits once p z ^\(z = 1) = p{z 
1 does not depend on V and p 2 ~i(l + ) = — b 2 /c 2 = p(l + ) 
and p 2 ^i(l~) = — 61/ci = p(l~) for V —> 00. 

In Fig. [1] we show the exact p vs. /3 p for different 
volumes V = L x L and r = 1/2. We compare such 
curves with Eq. ([8]) and also with Eq. ([2]) for the c 2 , 
61, and b 2 as above and a 2 //3 given by the best numeri- 
cal fitting to Eq. ([5|) (see Table I). For all the volumes, 
we observe a very good agreement between the exact p 
and the approximations, specially in the case where a 2 /(3 
is fitted. This latter demonstrates that the general Eq. 
(|2"|) - with the coefficients obtained from numerical sim- 
ulations - is an accurate procedure to calculate FOPTs 
order parameters. 

As a final analysis, we recall that the probability to 
be in phase x (for x = or r) is P x — w x /Z, with Z = 
wq + w r and w x the proper weight of phase x. Also, in 
the thermodynamic limit the term rz rV /(l + z' rV ) in the 
exact p, Eq. (O, is just the Heaviside function <d(z — 1) 
times the parameter r. Hence, it is the term z/(l + z) 
in Eq. §5§ that actually gives the density variation with 
respect to the fugacity. Thus, considering Eqs. (jU)-© 
close to the transition point - but still at the phase x - 
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TABLE I. The values of o 2 //3 used in Figure 1. 



L, where V — L x L 


a 2 //3 = rV 


02 //9 (numerical fitting) 


L = 4 


8 


10.016233 


L = 6 


18 


20.073368 


L = 8 


32 


34.164822 


L = 12 


72 


75.295494 



the probability and density, as function of z, read 



Pr 



1 



■exp[-pVf x 



1 + z 



rV : 



1 



(9) 



By isolating z in Eq. ([9]) one obtains P x (px), which 
gives (around z f=s 1) the probability to be in phase a; 
with the density value p x . Plots for two distinct V's and 
r = 1/2 are shown in the right inset of Fig. 1. It is inter- 
esting to observe that although P x (p x ) becomes broader 
for lower V's, even for a so small volume of V = 4 x 4 the 
density distributions of the two phases do not overlap if 
r = 1/2. On the other hand, by decreasing r, Pq and P r 
start to intersect each other. But r measures the jump 
in the order parameter, consequently how strong is the 
phase transition. Thus, a weaker phase transition leads 
to a larger overlap between the probability distributions 
of the order parameter close to the transition point. As 
it will be discussed in Sec. V, this result illustrates a 
general and important fact to set the method validity. 



IV. NUMERICAL APPLICATIONS 
A. Potts model 

Widely studied in statistical mechanics, both analyt- 
ically jrm20j and numerically HH^, the Potts model 
is a fine test for the present method, specially given its 
transition points are exactly known [8j. 

Consider each site of a regular lattice associated to 
a spin variable a, which assume the values 0, 1, . . . , q — 
1. If two adjacent sites have different (same) spin, their 
interaction energy is null (—J)- Thus, the Hamiltonian 
reads 



n = -jJ2 s ° 



(10) 



In the limit of very low temperatures, the system is con- 
strained to an ordered phase, becoming disordered as T 
increases. For any q, the order-disorder transition takes 
place at T c = l/ln[l + that in two dimensions is 
second-order for q < 4 and first-order for q > 5. An 
appropriate order parameter is 



q(Vmax/V) - 1 

Q-l 



(11) 



where V max is the volume occupied by the spins in a state 
a of lar gest po pulation and V — L 2 (in 2D) is the total 
volume 043. 




FIG. 2. The Potts model order parameter (j> versus T for 
(a) q — 20 and (b) q — 30, and different system sizes L. 
Continuous lines are from Eq. ||2j. The curves cross at To = 
0.5883(3) in (a) and at T = 0.5352(1) in (b). The insets show 
Tl versus 1/L 2 , where Tl is the temperature at the peak of 
X = -(d/dT)<j>. 



For the numerics we set q = 20 and q = 30, values 
which characterize strong first-order phase transitions. In 
Fig. [5]we show the order parameter, Eq. dill) , as function 
of T. We clearly see that in all cases <j> is well described by 
Eq. ([2]) (through proper parameters fitting). Moreover, 
in Fig. [5] the crossing points arc at Tq = 0.5883(3) for 
q = 20 and at T = 0.5352(1) for q = 30. Such values 
are corroborated by finite size scaling obtained from Tl 
versus 1/L 2 , with Tl the temperature at the peaks of 
the response function x — ~{9/dT)4> [HI, H2- Indeed, 
extrapolating the plots of Tl versus 1/L 2 (insets of Fig. 
IH, we find T = 0.5881(1) for q = 20 and T = 0.5353(3) 
for q — 30. The estimations are in excellent agreement 
with the exact values 1/ ln[l + V20] = 0.588349(. . .) and 
l/ln[l + V30] = 0.535248(...). 

As mentioned in the Sec. II. A, the mean energy and 
entropy, u and s, per site are also described by an ex- 
pression in the form of Eq. @. In Fig. [3] (Fig. 0]) we 
display the results for q = 20 (q = 30). Since entropy 
is not directly computed from standard thermodynamic 
averages, here we consider an indirect procedure, based 
on the transfer matrix method (for details refer, e.g., to 
4, 24]), to make the simulations and fit the coefficients 
of Eq. ^ in the case of s. We note that once more Eq. 
© indeed does describe quite well the relevant thermo- 
dynamic quantities and the crossing point (for a similar 
analysis for u and s, but considering a different approach, 
see Ref.[23j). We also can perform finite size scaling from 
the specific heat cy = (d/dT)u, plotting Tl versus 1/L 2 , 
for the Tl's the peak positions of the curves cyiT) for 
distinct Us QJ, [22[ (insets of Figs. [3] and [p. From the 
it's crossing and the Tl extrapolation we find, respec- 
tively, T = 0.5882(2) and T = 0.5882(1) for q = 20 
and T = 0.5353(4) and T = 0.5351(1) for q = 30, again 
consistent with the previous results. 

Finally, in Fig. [5] we display histograms of (j> for three 
different values of q = 10, 20, 30 at the coexistence. It 
illustrates that stronger the phase transition (i.e., higher 
the q's), lesser the overlap between the peaks (centered 
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FIG. 3. The Potts model with q = 20 and different L's. 
The continuous lines are from Eq. ([2]). (a) u versus T, with 
the crossing at T = 0.5882(2). The inset displays Tl versus 
l/L 2 , for Tl the temperature at the peak of the specific heat 
cv = (d/dT)u. (b) s versus T, with the inset showing the 
good coincidence of the crossing point. 
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-1.668 -1.662 -1.656 -1.65 -1.644 1.6558 1.6559 1.656 

FIG. 6. For the BL model, (a) p versus p for T — 0.3 and dis- 
tinct L's around the gas-LDL transition. All the continuous 
lines are properly obtained from Eq. ([2]). (b) The values of 
p — PL (for which (d/dp)p is maximum) versus l/L 2 . (c) A 
blow up of (a) around the crossing point and the respective 
small error bars. 
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T T 

FIG. 4. The same as Fig. 3, but for q = 30. For u, the curves 
cross at T = 0.5353(4). 




FIG. 5. For different Potts model q values and L = 10, his- 
tograms of the order parameter probability versus <j) at the 
coexistence. 



at the distinct phases cj> values) of the order parameter 
bimodal distribution. In particular, observe a larger over- 
lap for q = 10, for which we find from the present method 
a transition at T = 0.7016(9) (details not shown). The 
exact value is l/ln[l + V^0] = 0.701231(. . .). Therefore, 
although still good, it is not so accurate as the previous 
examples. 



B. Bell-Lavis model 

The Bell-Lavis (BL) @] is a lattice gas model able to 
reproduce liquid polimorphism and water-like anomalies. 
It is defined on a triangular lattice where each site is 
characterized by its occupation (er) and orientation (r) 
states. Whether the site i is or is not occupied by a water 
molecule, cr, = 1 or crj = 0, respectively. Furthermore, 
if the site i has an "arm" which is (is not) inert towards 
the adjacent site j, then — {t? = 1). 

Two nearest neighbor molecules i and j interact via 
a van-der-Waals energy —e v dw They also form a hy- 
drogen bond (of energy — eub) when t^Tj 1 = 1. So, in 
the grand-canonical ensemble the BL is described by the 
Hamiltonian 

H= - ^ <J t (J j ( e hb T* 3 + e vdw ) - n^ai. (12) 

<i,j> i 

The van-der-Waals interaction favors an increasing in the 
lattice density (a proper order parameter), whereas the 
hydrogen bond tends to form sublattices for which the 
molecules have opposite orientations. 

For C = e v dw/^hb < 1/3, the system presents three 
stable phases, named: gas, low-density- liquid (LDL), and 
high-density-liquid (HDL). For low n, the system is in 
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the gas phase. By increasing p we go through the gas- 
LDL and then through the LDL-HDL phase transitions. 
At zero temperature both gas-LDL and LDL-HDL are 
of first order, taking place at p c = — 3 (1 + Q/2 and 
p c = —6 £, respectively. For T > 0, the gas-LDL remains 
first-order (ending up in a tricritical point if C = 1/10), 
but the LDL-HDL becomes second-order 0, Oj, [H, H . 

For T — 0.3, in Fig. [B] (a) we plot p versus p around 
the gas-LDL coexistence. As previously, the isotherms 
are well described by Eq. @, with a crossing occurring 
at p = —1.6559(1) for p ss 0.516(2). Such p is close to 
1/2, the exact result at T = (understood recalling that 
at the coexistence both gas (p = 0) and LDL (p = 2/3) 
phases have equal weight and the LDL has degeneracy 
O.LD l = 3). In Fig. [6] (b) we display pl (the p for which 
(d/dp)p is a maximum) versus 1/L 2 . By taking the ther- 
modynamic limit L — ¥ oo, we find the extrapolated value 
of po = —1.6560(1), in excellent agreement with the es- 
timate in Fig. [6] (a). Finally, Fig. [6] (c) is a blow up of 
Fig. [5] (a) in the vicinity of the phase coexistence. The 
observed small error bars illustrate the good accuracy of 
the present method in locating the transition point. 

We should stress that although this is the only instance 
where we present a more detailed error analysis, in all the 
other examples the error bars are likewise small. 



C. Associating lattice-gas (ALG) model 

Similarly to the BL, the symmetric associating lattice- 
gas (ALG) model [l(| EH can display liquid polimorphism 
and water anomalies. It is also defined on a triangular 
lattice, where each site is described by an occupation 
(cr) and orientation (r) state. But an important differ- 
ence from the BL is that an energetic punishment exists 
when a hydrogen bond is not formed. Two first neighbor 
molecules have an interaction energy of — v {—v + 2u) if 
there is (there is not) a hydrogen bond between them. 
The Hamiltonian therefore reads 

H = 2u £ ^[(l-^u))-^]-^^. (13) 

<i,j> i 

The ALG presents a gas and two liquid, LDL and HDL, 
phases. In particular, for the LDL phase 3/4 of the lattice 
is filled by water molecules forming the maximum num- 
ber of hydrogen bonds [11| . Another relevant distinction 
from the BL model is that here both gas-LDL and LDL- 
HDL transitions remain first-order for T ^ 0. At T = 0, 
the discontinuous transitions occur at p/v = —2 (gas- 
LDL) and p/v = -6 + 8u/v (LDL-HDL). 

Around the transitions gas-LDL, Fig. [3(a), and LDL- 
HDL, Fig. |S] (a), we plot the order parameter versus p 
for T = 0.20 and different Us. For the former, we simple 
take p as the order parameter. However, since for LDL- 
HDL the density is never null, we set as the order param- 
eter <j> = (4/9—3). Both cases are completely described by 
Eq. @, with the crossing occurring at po = —2.0000(2), 




-2.04 -2.02 -2 -1.98 -1.96 0.25 - 0.5 0.75 



FIG. 7. For the ALG model and the gas-LDL transition, (a) 
the density p versus p for distinct volumes and T — 0.2, (b) 
the linear scaling of pl versus 1/L 2 , and (c) the probability 
density P p versus p for L = 12 at the coexistence. Results for 
V — 4 x 2 are exact. 
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M- <j) 

FIG. 8. Similar to Fig. but for the LDL-HDL transition 
and the order parameter (j> (see main text). 



gas-LDL, and p = 2.0000(2), LDL-HDL. These esti- 
mates (within the numerical uncertainties) are identical 
to their exact values at T = 0, a particularity of the 
ALG model. However, by increasing more the tempera- 
ture, the pqs start to change as well. For instance, for 
gas-LDL (which has a shorter coexistence line than that 
for LDL-HDL [11]) p = -1.9986(2) at T = 0.3 @|. On 
the other hand, for LDL-HDL it is necessary T > 0.5 for 
a sensible departure of po from its value of 2 at T = 
(see also Sec. V), e.g., for T = 0.6, we have found that 
po = 1-9970(5) (results not shown). 

In Figs. [7] and [S] (b) we plot pl versus 1/L 2 , for 
Pl the p for which the respective response functions, 
{d/dp)p and (d/dp)4>, present a maximum. Extrapo- 
lating to the thermodynamic limit we get the estimates 
po = -1.9999(1) and = 2.0000(1), values very close 
to those from the crossing calculation. Lastly, the bi- 
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FIG. 9. (a) For the 3NN hard-core gas model, a filled site 
(•) hinders neighbor sites (x) to be occupied, (b) This cor- 
responds to an spatial exclusion region interaction (dashed 
circle), (c) The regular structure in the maximum possible 
filling configuration has (d) an unitary cell with four vacant 
sites (in a total of five), resulting in a density of p ma x = 1/5. 



modal density probability distributions, P p and P^, are 
shown Figs. [7] and [5] (c). They present a very flat val- 
ley between the peaks (with each peak associated to an 
individual phase at the coexistence). 



D. Hard-core gas lattice model 

As a last example, let us assume a 2D square lattice for 
the hard-co re g as model introduced in |27[ and recently 
revisited in [13 • The interaction is entirely given by 
an exclusion range: a particle in a certain site prevents 
the occupation of all the surrounding sites. In the so 
called 3NN version, the one analyzed here, the excluded 
sites are those shown in Fig. [§] (a)-(b). In this case, the 
lattice maximum possible filling is displayed in Fig. [9] 
(c), resulting in a density of p max = 1/5 (Fig. [§] (d)). 

The sole control parameter is an effective chemical po- 
tential ft, which determines the total number of particles 
in the system. Hence, temperature is not defined, char- 
acterizing an athermal problem. For a fixed ft, the proba- 
bility to have n particles in a lattice of volume V is given 
by p(n) = a n exp[ftn]/Z, for Z = Y^Zo^ a n exp[ftn], 
Nmax = V/5, and a n the number of distinct configura- 
tions of n particles allowed by the hard-core potential. By 
decreasing (increasing) ft, the system density decreases 
(increases). For lower ft's, the particles are basically ran- 
domly distributed - obeying the above restrictions - con- 
stituting a fluid phase. For higher ft's, the system starts 
to present a certain ordering so to accommodate larger 
numbers of particles (up to a maximum of N max ). The 
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FIG. 10. (a) A basic cell (of 5 x 5 sites) for the 3NN hard-core 
gas model (here illustrated in the maximum filling condition). 
The full lattice is formed by juxtaposing basic cells, (b) The 
basic cell sites labeling, with a;+5fc = a; for i = 1, . . . , 5 and 
k = 0, . . . , 4. So, there are five distinct ways to name a basic 
cell, but once one is assumed, it should be used throughout. 

transition fluid-ordering is of first-order (actually, in the 
ordering regime there are two phases related to each other 
by a chiral transformation |2a - [30l |). 

Since the usual density p only vanishes for ft — > — oo 
and the phase transition takes place for a finite ft, a more 
appropriate order parameter (j> should be considered. For 
so, we follow the approach in [29| and take the full lattice 
as composed of basic cells (sublattices) of 5 x 5 sites each, 
Fig. [TO] (a). The sites in each row of a basic cell are 
labeled as a,i (i — 1,...,5) and in total there are five 
different ways (k = 0, . . . , 4) to name it (sec Fig. [TO] (b)). 
Thus, from such construction one can set </> as in [29( 
(just using a slight different notation), or 

0HI&=3-&=o|>, h = ^ £ <-<|. 

(14) 

Above, (...) denotes average over all the lattice basic 
cells. Also, for a chosen labeling k, n\. denotes the 
number of particles in the sites named of a basic 
cell. According to this definition, at the maximum filling 
fc=o = and 4>k=3 = 1, so = 1. On the other hand, at 
low densities (fluid phase) 4>k tends to zero, consequently 
it does 4>. 

The numerical simulation procedure is essentially that 
described in Sec. II. B. The only small differences are: 
(i) instead to define the replicas at distinct tempera- 
tures, they are defined at distinct /i's; (ii) the occupation 
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FIG. 11. For the 3NN hard-core gas model, the order param- 
eter <j) versus the effective chemical potential jl for different 
L's. Continuous lines are from Eq. The curves cross at 
flo — 3.6741(8). In the inset fiL versus 1/L 2 , where fiL is the 
position of the peak of the "susceptibility" x = (d/dji)(j). 



state of a site (observing the exclusion rule) is changed 
according to min{l, exp[±/2]}, where the signal + (— ) 
is taken if the site is initially empty (occupied); finally 
(iii) the exchange of configuration between two replicas, 
say A and B, is performed according to the probability 
min{l, exp[(p,B — Aa)AA]}, with AN the difference of 
the number of particles of A and B. 

We have simulated the model for system sizes rang- 
ing from L = 25 to L = 40, shown in Fig. [TTJ Note 
that all curves are very well described by Eq. © , whose 
crossing point occurs at (// o ,0o) = (3.6741(8), 0.835(6)). 
Thus, even for an uncommon (but appropriate) definition 
for the order parameter, Eq. (|14|) , it is properly repre- 
sented by our general function W. In the inset of Fig. [11] 
we plot fiL versus 1/L 2 , where /*l is the position of the 
peak of (d/dp)4>. In the thermodynamic limit we find 
the value p, — 3.6758(9), in agreement with the crossing 
estimate and with the values 3.6762(1) in [2_8J and 3.6746 
in 29] . In Fig. [T^] the probability distribution histogram 
of the order parameter for L = 30 and at the coexis- 
tence presents a low valley between the peaks of the two 
coexisting phases. 



V. THE METHOD APPLICABILITY 

So far we have discussed five distinct examples: a 
prototype thermodynamic system, three representative 
lattice-gas models, and an interesting athermal prob- 
lem. We have found that the present approach is able 
to describe FOPTs in all the different situations studied. 
Thus, a relevant issue is to enquire to what extent the 
method can give good results. 

To address it, we first recall that the approach key 
point is the actual form of Eq. (fT]), i.e., to write Z as a 
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FIG. 12. For the 3NN hard-core gas model and L = 30, 
histogram of the order parameter probability P$ at the coex- 
istence. 
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FIG. 13. (a) In an arbitrary phase space £ x f, example of 
a small region Q around a point of coexistence of N phases 
(here, M = 6). (b) Schematics of F versus £ for ( = (* . If 
^ > the system is in the phase n = 5 and the dashed 
line represents the functional form of F\ extended into such 
region. From the inclination of J5, the FOPT 1-5 would be 
stronger in (ii) than in (i). So, Ai, 5 (£') = Fi(£') - F 5 (f) 
would be larger in the case (ii). 

sum of exponentials, where each term is uniquely associ- 
ated to a particular phase. In other words, there are no 
terms involving overlapping between coexisting phases. 
The rigorous analysis in [lij HU show that this decom- 
position is generally valid at low temperatures (see also 
the Appendix). 

To understand in a more physical ground why of a such 
structure for the partition function, one might consider 
the following heuristic arguments (in the specific case we 
are close to a FOPT point): 

(i) Assume in phase space a small region Q encompass- 
ing a point of coexistence of A/" phases. Moreover, 
for any n — 1 , . . . , Af, let f2„ be a portion of f2 
corresponding to the phase n (e.g., Fig. [TBI (a)). 

(ii) Within any ft n , F n — — \n[Z]/(3 is the phase n free 
energy. Then, suppose (at least formally) that in 
f2„ we can write Z w exp[— (3F n ](l + Z n ), with a 
proper Z n being very small in such region (in fact, 
we must have \Z n /j3\ -C F n ). 

(iii) Thus, in each f2„: Z w exp[— j3F n ] + a small term. 
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Therefore, a tentative partition function for the 
whole SI can be 2 k S^Li ex P[ — PFn] (for sim- 
plicity neglecting possible degeneracies a n ). This 
is just Eq. (JTJ. 

(iv) But for the above to hold, a consistency condition is 
required. Notice that Z n = 2 n /^„ exp[-/3A„< in ], 
with A. n i >n = F n i — F n representing the difference 
between the free energy of phases n' and n. So, Z n 
to be small in Sl„ implies that Vn' ^ n, /3A„> )n is 
large in Sl„. 

From the above reasoning we reach the desired general 
validity condition for Eq. (TTJ (and so for W in Eq. ©), 
namely, 

/3A„/ jn ^> 1 in fl n for any n and for all n' 7^ n. (15) 

Note that it explains why Eq. (J]) is always valid at low 
T's (at least close to phase transitions). Indeed, in such 
case even if the A's are not large, the product (3 A can be 
very large if the temperature is sufficiently small. 

But Eq. (fT5|) is also true if the A's themselves are 
large (of course, with T not too high). As illustrated in 
Fig. [T3] (b), this is the case in strong FOPTs, i.e., for 
the system displaying large discontinuities in the slop of 
the F's across £* or, equivalently, for a large jump in the 
value of the order parameter (for instance, as determined 
by the quantity r in the example of Sec. III). 

Lastly, there is a practical and computationally inex- 
pensive test to check the above relations. Equation ([15)) 
implies in very high entropic barriers across the tran- 
sition point. Hence, even considering ubiquitous ther- 
modynamic fluctuations (e.g., for finite systems) around 
such point it would be much more probable the order 
parameter to assume values typical of the single phases 
(corresponding to the F's minima) than to present val- 
ues in between (implying the system to cross the high F's 
regions). So, one could calculate the probability distri- 
bution histogram for the order parameter at the coexis- 
tence condition, for which the peaks relate to the distinct 
phases. The verification of Eq. (|T5|) would result in well 
separated peaks, not overlapping each other. Indeed, this 
is exactly the case in the examples here, as observed in 
Figs, m (right inset), [51 0(c), M (c), and[TJ 

As a final illustration, we come back to the ALG model 
of Sec. III.C. However, we set T = 0.5, which is 2.5 times 
higher than the value in Figs. [7] and [8] This case is inter- 
esting because now the gas-LDL (LDL-HDL) transition, 
Fig. [TJ](Fig. US]), is well described by Eq. @ only closer 
to the phases coexistence point. Indeed, compare the fit- 
ting quality and the p range considered in Figs. [7J and 
IT4l and in Figs. [5] and [15] Important to emphasize that 
unlike Fig. [JJ(c), in Fig. [U](c) the probability distribu- 
tion is not really a flat valley between the peaks for the 
gas-LDL. On the other hand, similar to Fig. [5] (c), for 
LDL-HDL the two peaks in Fig. [T3] (c) do not intersect. 

Nevertheless, the range in which W is valid is still large 
enough to allow a proper characterization of the FOPT 




Li p 



FIG. 14. For the ALG model gas-LDL transition at T = 0.5, 
(a) the density p versus p for distinct L's. Continuous lines 
are obtained from Eq. @. (b) The scaling plot of pl (for 
which x = ((p 2 ) — {p) 2 W i s a maximum) versus 1/L 2 . (c) 
The probability density of p at the coexistence for L = 12. 




H if 



FIG. 15. The same as in Fig. but for the LDL-HDL 
transition. Here p,L is the peak position of \ — (p/dfj,)ij>. 



in the gas-LDL case (for LDL-HDL, the range is even 
larger, see Fig. [15]). From the crossing curves in Fig. 
Il"4l we get the estimate ^0 = —1-9360(5), in agreement 
with fiQ = —1.9365(5) from the position /i^ of the peak 
of x = ((/o 2 ) - (p) 2 )V, Fig. [JJ (b). We observe that 
since the fitting of 4> is not so good for a broader interval 
of jti's, we prefer to calculate \ = ((p 2 ) — (p) 2 )^ (thus 
numerically more reliable) instead to define \ — {9/dfi)(p. 
For the LDL-HDL case, the curves crossing leads to the 
estimate [iq = 1.9995(5), very close to the estimate /io = 
1.9990(5), obtained from the peak of \ = {d/dp)<j). 



VI. REMARKS AND CONCLUSION 

In this contribution we have clarified the main mathe- 
matical aspects, discussed the applicability condition and 
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extended the instances of usage of a recent proposed [|[ 
approach to treat FOPTs, which can be summarized as 
the following. From a special decomposition for the par- 
tition function Eq. ((T]) - valid close to a FOPT when- 
ever Eq. (|T5|) holds - one can derive an analytical ex- 
pression W, Eq. ((2j, which depends on few coefficients. 
By using simple numerical simulations to determine these 
coefficients, W is able to describe quite well relevant ther- 
modynamic quantities, like order parameter, energy, en- 
tropy, etc, around the transition point. In addition, there 
is a point where all curves W (irrespect of L) cross. By 
considering relatively small system sizes L, the crossing 
can give the transition point thermodynamic value. 

As it should be, the procedure agrees with other effi- 
cient schemes available [1, [IH, [r| [29|. However, it 
has the advantage of being general, inexpensive from the 
computational point of view and to yield response func- 
tions x ( e -g-: compressibility and specific heat) in a rather 
direct way (analytically). 

The method validity condition can be tested from plots 
of multimodal probability distributions of the order pa- 
rameter at the coexistence, calculated from straightfor- 
ward simulations. A non-overlapping of the peaks (as- 
sociated to the individual phases) indicates that it can 
be satisfactorily applied. In fact, such condition extends 
the protocol, originally derived pj for the situation of low 
temperatures. 

The method has been tested and shown to work fine 
for several lattice problems, including an exact solvable 
and the relevant Potts, Bell-Lavis, ALG, and athermal 
hard-core gas, models. But certainly, a natural question 
is if FOPTs taking place in off-lattice systems - with 
much larger phase spaces - could be studied in a simi- 
lar fashion. In this respect, we first observe that some 
aspects of continuous systems displa ying FOPTs at low 
T's, (e.g., as for polymers in Ref. [32(, analyzable in 
terms of finite-size scaling (33|) can be described by lat- 
tice models [34]]. Obviously, in such cases the method 
could be directly applied. Second, note that the decom- 
position for the partition function Z, when valid, does not 
make restrictions regarding lattice or off-lattice systems. 
Thus, the only issue would be the use (in a continuous 
phase space) of a proper simulation sampling procedure 
to fit the parameters in Eq. For instance, there are 
some implementations of the PT for off-lattices systems 
(see, e.g., Refs. (35|). With such implementations the 
approach should hold in the same way. 

The study of off-lattice [35| and polymers systems dis- 
playing FOPTs [36[ is presently an ongoing work and will 
be reported in the due course. 
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FIG. 16. Examples of possible phase diagrams. By varying 
the control parameter £ in an interval y = £ — £*, it is as- 
sumed that other intensive quantities, f, are properly fixed, 
(a) Along the displayed separation line between two phases, 
any £ always will allow a £*. (b) On the other hand, for the 
shown triple point, £ should be set to £*. 



Appendix A: The derivation of W and some of its 
properties 

Consider a finite system at a low temperature and pre- 
senting N coexisting stable phases (the meaning of "low" 
here is discussed in Section V). It has been rigorously 
shown [jjl, Hl| that the problem partition function is well 
described by (with /3 = (kT)^ 1 and V the volume) 



z = zZ 

n=l 



a n exp[-/3Vf n ] + Z„ 



(Al) 



Z U nst is associated to the possible existence of unstable 
phases - but which are exponentially damped, so negligi- 
ble in Eq. (|A1|) - and /„ is the n-th phase [n = 1, . . . , N) 
free energy per volume The degeneracy parameters 
(or weights) a's result from eventual symmetries, so that 
a n > 1 would be due to distinct spatial configurations 
leading to a same phase n. 

Now, let £ to be a proper phase transition control pa- 
rameter, which we shall vary. From Eq. (jAll) we define 



W 



with 



1 d 



HZ] = 



9r> 



2^71=1 



On g n exp[-/3Vf„ 



22 n =i a » ex Pr 



Wfn] 



T 



(A2) 



(A3) 



Note this the definition of W is usually the start point 
to calculate relevant order parameters. 

The system may have many other intensive parame- 
ters, which we generally denote by £. So, suppose £ kept 
fixed at proper values, such that the Af phases coexis- 
tence takes place at £ = £* (usually with £* depending 
on C, see Fig. [T6|) . In this case, the f n 's are single func- 
tions of £ and for y = £ — £* w we consider a first 
order series expansion (possible due to the existence of 
smooth representations for the / n 's [3l|): f n ~ f*+fn Vi 
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with / n (£*) = /* V n (because the coexistence) and 
/£ = (0/ n /00| c=c .. It leads to 



f" 

■I 1 1 



/* dT 



(A4) 



where, obviously, dT/d£ — 1 for £ = T and zero other- 
wise. In deriving Eq. (|A4[) . if the control parameter £ 
is the temperature, one also must assume 1/T w 1/T*, 
a reasonable approximation provided VF is calculated for 
|T — T*\ small (for how small in practice, see the numer- 
ical examples along the paper). 

Next, for a n = {f* - f[*) V$ [or a n = (f n * - /(*) V(3* 
if f = T], b n = (a n /ai)g n , and c„ = (a n /ai), we get 



& i + J2n=2 b n exp[-a„y] 



W 



1 + YH=2 c « exp[-a„y] 



(A5) 



Above, the coefficients a n , b n and c„ are independent 
on the control parameter and only the c n 's depend (lin- 
early) on the volume. Hence, at the coexistence (y = 0) 
W 7^ W(V) and all the curves W versus £, regardless 
of V, must cross at £ = £*. Thus, Eq. (|A5I) not only 
describes generic thermodynamic quantities, but also 
yields the thermodynamic limit estimate for the tran- 
sition point. Furthermore, by taking derivatives of Eq. 
(IA5|) . one can obtain response functions and susceptibil- 
ities. 

Finally, at £ = £* either from Eq. (|A2j) or from Eq. 
~TB W reads 



f* r)T 

71—1 ^ 



(A6) 



with 



Pn = 



TFT 



(A7) 



Moreover, suppose the /'s ordered such that f* = /{* 



f * < f 



m+1 



<--.<f' k -i </i* = ••• = /£ = A*. 

Assuming y small, so we can consider Eq. (|A5I) . if we 
take V —¥ oo with y positive (case +) or y negative (case 
— ), for v+ = 1, u+ = rruv- = k,U- = J\f we find (note a 
little misprint in Ref. |5() 



W+ = 



t—tr, 



bn 



Eri—u±- 
n—v± 



= f „ _ r dT 

T* d£ ' 



(A8) 



Equations (|A6|) - (|A8|) give the discontinuity of the order 
parameter across the phase transition in the thermody- 
namic limit. Numerically, such discontinuity is obtained 
from the coefficients 6's and c's once, from Eq. (IA8|) . we 
can write W + — b\jc\ and W_ — bj^/cj^f. In particular, 
for k = Af (m = 1), W+ (W_) is given in terms of the sole 
phase which is immediately to the right (left) of £*. Also, 
the number of equal a„'s correspond, at least in a first 
order approximation, to the number of phases which co- 
exist over the line £ in the vicinity of £* . This fact can be 
used to locate coexisting phase lines and triple points (an 
application for the present method to appear elsewhere). 
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